coverage.pool <- colMeans(ifelse(lllist.pool[[index]] > b1 | ullist.pool[[index]] < b1, 0, 1))
plot(truecor, coverage.fe, pch=19, type="l", col = "#80b1d3", lwd = 3, ylim=c(0, 1), xlim=c(0, 1),
xlab=expression(rho), ylab="Coverage",
main=paste(unitloop, "units, ", Nperloop, "observations per unit"))
lines(truecor, coverage.re, pch=19, col = "#fb8072", lwd = 3, lty=2)
lines(truecor, coverage.pool, pch=19, col = "#bebada", lwd = 3, lty=3)
index <- index + 1
}
}
dev.off()
# Plot coverage (how often do our 95% CI's reject beta= truth?):
## FE=solid line, RE=dashed line, pooled=dotted line
#quartz(width=9, height=8)
pdf("CL_Coverage_Sluggish.pdf", width=9, height=8)
par(mfrow=c(length(units), length(Nperunit)), mar=c(4, 4, 2, 2), las=1)
index <- 1
for (unitloop in units) {
for (Nperloop in Nperunit) {
coverage.fe <- colMeans(ifelse(lllist.fe[[index]] > b1 | ullist.fe[[index]] < b1, 0, 1))
coverage.re <- colMeans(ifelse(lllist.re[[index]] > b1 | ullist.re[[index]] < b1, 0, 1))
coverage.pool <- colMeans(ifelse(lllist.pool[[index]] > b1 | ullist.pool[[index]] < b1, 0, 1))
plot(truecor, coverage.fe, pch=19, type="l", col = "#80b1d3", lwd = 3, ylim=c(0, 1), xlim=c(0, 1),
xlab=expression(rho), ylab="Coverage",
main=paste(unitloop, "units, ", Nperloop, "observations per unit"))
lines(truecor, coverage.re, pch=19, col = "#fb8072", lwd = 3, lty=2)
lines(truecor, coverage.pool, pch=19, col = "#bebada", lwd = 3, lty=3)
index <- index + 1
}
}
dev.off()
pdf("CL_RMSE_Sluggish.pdf", width=9, height=8)
par(mfrow=c(length(units), length(Nperunit)), mar=c(4, 4, 2, 0.5), las=1)
index <- 1
for (unitloop in units) {
for (Nperloop in Nperunit) {
rmse.fe <- sqrt(colMeans((coeflist.fe[[index]] - b1)^2))
rmse.re <- sqrt(colMeans((coeflist.re[[index]] - b1)^2))
rmse.pool <- sqrt(colMeans((coeflist.pool[[index]] - b1)^2))
plot(truecor, rmse.fe, pch=19, type="l", col = "#80b1d3", lwd = 3, ylim=c(0, 1), xlim=c(0, 1),
xlab=expression(rho), ylab="RMSE",
main=paste(unitloop, "units, ", Nperloop, "observations per unit"))
lines(truecor, rmse.re, col = "#fb8072", lwd = 3, lty=2, pch=19  )
lines(truecor, rmse.pool, col = "#bebada", lwd = 3, pch=19, lty=3)
index <- index + 1
}
}
dev.off()
#quartz.save("CL_RMSE_Sluggish", type = "pdf")
# Plot mean bias:
## FE=solid line, RE=dashed line, pooled=dotted line
#quartz(width=9, height=8)
pdf("CL_Bias_Sluggish.pdf", width=9, height=8)
par(mfrow=c(length(units), length(Nperunit)), mar=c(4, 4, 2, 0.5), las=1)
index <- 1
for (unitloop in units) {
for (Nperloop in Nperunit) {
mab.fe <- colMeans(coeflist.fe[[index]] - b1)
mab.re <- colMeans(coeflist.re[[index]] - b1)
mab.pool <- colMeans(coeflist.pool[[index]] - b1)
plot(truecor, mab.fe, pch=19, type="l", col = "#80b1d3", lwd = 3, ylim=c(0, 1), xlim=c(0, 1),
xlab=expression(rho), ylab="Bias",
main=paste(unitloop, "units, ", Nperloop, "observations per unit"))
lines(truecor, mab.re, pch=19, col = "#fb8072", lwd = 3, lty=2)
lines(truecor, mab.pool, pch=19, col = "#bebada", lwd = 3, lty=3)
index <- index + 1
}
}
dev.off()
#quartz.save("CL_Bias_Sluggish", type = "pdf")
# Plot standard deviation:
## FE=solid line, RE=dashed line, pooled=dotted line
#quartz(width=9, height=8)
pdf("CL_SD_Sluggish.pdf", width=9, height=8)
par(mfrow=c(length(units), length(Nperunit)), mar=c(4, 4, 2, 0.5), las=1)
index <- 1
for (unitloop in units) {
for (Nperloop in Nperunit) {
sd.fe <- colSds(coeflist.fe[[index]])
sd.re <- colSds(coeflist.re[[index]])
sd.pool <- colSds(coeflist.pool[[index]])
plot(truecor, sd.fe, pch=19, type="l", col = "#80b1d3", lwd = 3, ylim=c(0, 1), xlim=c(0, 1),
xlab=expression(rho) ylab="SD",
main=paste(unitloop, "units, ", Nperloop, "observations per unit"))
lines(truecor, sd.re, pch=19, col = "#fb8072", lwd = 3, lty=2)
lines(truecor, sd.pool, pch=19, col = "#bebada", lwd = 3, lty=3)
index <- index + 1
}
}
dev.off()
#quartz.save("CL_SD_Sluggish", type = "pdf")
# Plot power (1-TypeII error):
## FE=solid line, RE=dashed line, pooled=dotted line
#quartz(width=9, height=8)
pdf("CL_Power_Sluggish.pdf", width=9, height=8)
par(mfrow=c(length(units), length(Nperunit)), mar=c(4, 4, 2, 0.5), las=1)
index <- 1
for (unitloop in units) {
for (Nperloop in Nperunit) {
power.fe <- colMeans(ifelse(lllist.fe[[index]] > 0 | ullist.fe[[index]] < 0, 1, 0))
power.re <- colMeans(ifelse(lllist.re[[index]] > 0 | ullist.re[[index]] < 0, 1, 0))
power.pool <- colMeans(ifelse(lllist.pool[[index]] > 0 | ullist.pool[[index]] < 0, 1, 0))
plot(truecor, power.fe, pch=19, type="l", col = "#80b1d3", lwd = 3, ylim=c(0, 1), xlim=c(0, 1),
xlab=expression(rho), ylab="Power",
main=paste(unitloop, "units, ", Nperloop, "observations per unit"))
lines(truecor, power.re, pch=19,  col = "#fb8072", lwd = 3, lty=2)
lines(truecor, power.pool, pch=19, col = "#bebada", lwd = 3, lty=3)
index <- index + 1
}
}
dev.off()
#quartz.save("CL_Power_Sluggish", type = "pdf")
# Plot coverage (how often do our 95% CI's reject beta= truth?):
## FE=solid line, RE=dashed line, pooled=dotted line
#quartz(width=9, height=8)
pdf("CL_Coverage_Sluggish.pdf", width=9, height=8)
par(mfrow=c(length(units), length(Nperunit)), mar=c(4, 4, 2, 0.5), las=1)
index <- 1
for (unitloop in units) {
for (Nperloop in Nperunit) {
coverage.fe <- colMeans(ifelse(lllist.fe[[index]] > b1 | ullist.fe[[index]] < b1, 0, 1))
coverage.re <- colMeans(ifelse(lllist.re[[index]] > b1 | ullist.re[[index]] < b1, 0, 1))
coverage.pool <- colMeans(ifelse(lllist.pool[[index]] > b1 | ullist.pool[[index]] < b1, 0, 1))
plot(truecor, coverage.fe, pch=19, type="l", col = "#80b1d3", lwd = 3, ylim=c(0, 1), xlim=c(0, 1),
xlab=expression(rho), ylab="Coverage",
main=paste(unitloop, "units, ", Nperloop, "observations per unit"))
lines(truecor, coverage.re, pch=19, col = "#fb8072", lwd = 3, lty=2)
lines(truecor, coverage.pool, pch=19, col = "#bebada", lwd = 3, lty=3)
index <- index + 1
}
}
dev.off()
#quartz.save("CL_Coverage_Sluggish", type = "pdf")
# Plot Optimism
## FE=solid line, RE=dashed line, pooled=dotted line
#quartz(width=9, height=8)
pdf("CL_Optimism_Sluggish.pdf", width=9, height=8)
par(mfrow=c(length(units), length(Nperunit)), mar=c(4, 4, 2, 0.5), las=1)
index <- 1
for (unitloop in units) {
for (Nperloop in Nperunit) {
sd.fe <- colSds(coeflist.fe[[index]])
sd.re <- colSds(coeflist.re[[index]])
sd.pool <- colSds(coeflist.pool[[index]])
se.fe <- colMeans(selist.fe[[index]])
se.re <- colMeans(selist.re[[index]])
se.pool <- colMeans(selist.pool[[index]])
opt.fe <- se.fe/sd.fe
opt.re <- se.re/sd.re
opt.pool <- se.pool/sd.pool
plot(truecor, opt.fe, pch=19, type="l", col = "#80b1d3", lwd = 3, ylim=c(0, 1), xlim=c(0, 1),
xlab=expression(rho), ylab="Overconfidence",
main=paste(unitloop, "units, ", Nperloop, "observations per unit"))
lines(truecor, opt.re, pch=19, col = "#fb8072", lwd = 3, lty=2)
lines(truecor, opt.pool, pch=19, col = "#bebada", lwd = 3, lty=3)
index <- index + 1
}
}
dev.off()
#quartz.save("CL_Optimism_Sluggish", type = "pdf")
# end of file.
library(lme4)
library(MASS)
library(matrixStats)
## FE=solid line, RE=dashed line, pooled=dotted line
#quartz(width=9, height=8)
pdf("CL_Optimism_Sluggish.pdf", width=9, height=8)
par(mfrow=c(length(units), length(Nperunit)), mar=c(4, 4, 2, 0.5), las=1)
index <- 1
for (unitloop in units) {
for (Nperloop in Nperunit) {
sd.fe <- colSds(coeflist.fe[[index]])
sd.re <- colSds(coeflist.re[[index]])
sd.pool <- colSds(coeflist.pool[[index]])
se.fe <- colMeans(selist.fe[[index]])
se.re <- colMeans(selist.re[[index]])
se.pool <- colMeans(selist.pool[[index]])
opt.fe <- se.fe/sd.fe
opt.re <- se.re/sd.re
opt.pool <- se.pool/sd.pool
plot(truecor, opt.fe, pch=19, type="l", col = "#80b1d3", lwd = 3, ylim=c(0, 1), xlim=c(0, 1),
xlab=expression(rho), ylab="Overconfidence",
main=paste(unitloop, "units, ", Nperloop, "observations per unit"))
lines(truecor, opt.re, pch=19, col = "#fb8072", lwd = 3, lty=2)
lines(truecor, opt.pool, pch=19, col = "#bebada", lwd = 3, lty=3)
index <- index + 1
}
}
dev.off()
pdf("CL_SD_Sluggish.pdf", width=9, height=8)
par(mfrow=c(length(units), length(Nperunit)), mar=c(4, 4, 2, 0.5), las=1)
index <- 1
for (unitloop in units) {
for (Nperloop in Nperunit) {
sd.fe <- colSds(coeflist.fe[[index]])
sd.re <- colSds(coeflist.re[[index]])
sd.pool <- colSds(coeflist.pool[[index]])
plot(truecor, sd.fe, pch=19, type="l", col = "#80b1d3", lwd = 3, ylim=c(0, 1), xlim=c(0, 1),
xlab=expression(rho), ylab="SD",
main=paste(unitloop, "units, ", Nperloop, "observations per unit"))
lines(truecor, sd.re, pch=19, col = "#fb8072", lwd = 3, lty=2)
lines(truecor, sd.pool, pch=19, col = "#bebada", lwd = 3, lty=3)
index <- index + 1
}
}
dev.off()
colMeans(ifelse(lllist.re[[3]] > b1 | ullist.re[[3]] < b1, 0, 1))
lllist.re[[3]]
View(lllist.re)
test <- lllist.re[[3]]
View(test)
summary(test[11])
missing_ll <- is.na(lllist.re[[3]])
View(missing_ll)
tabulate(missing_ll)
missing_ll <- as.numeric(is.na(lllist.re[[3]]))
tabulate(missing_ll)
count(missing_ll)
missing_ll
sum(missing_ll)
test <- ullist.re[[3]]
missing_ul <- as.numeric(is.na(ullist.re[[3]]))
sum(missing_ul)
View(test)
colMeans(ifelse(lllist.re[[3]] > 0 | ullist.re[[3]] < 0, 1, 0))
colMeans(ifelse(lllist.re[[3]] > b1 | ullist.re[[3]] < b1, 0, 1), na.rm = T)
View(test)
View(coeflist.re)
View(coeflist.re)
View(coefmat.fe)
View(coefmat.re)
View(coefmat.re)
View(semat.re)
View(coeflist.pool)
View(coeflist.pool)
test <- coeflist.re[[3]]
View(test)
test <- selist.re[[3]] ## no missing values
View(test)
load("C:/Users/alikagalwala/Dropbox/CSTS-Panel Papers/JOP RR/Replication Files/Clark and Linzer 2015/Sluggish/2000sims_sluggish.RData")
# Clark and Linzer replication with additional statistics from the MC results
# ----------------------------#
set.seed(1138183)
library(lme4)
library(MASS)
library(matrixStats)
pdf("CL_Optimism_Sluggish.pdf", width=9, height=8)
par(mfrow=c(length(units), length(Nperunit)), mar=c(4, 4, 2, 0.5), las=1)
index <- 1
for (unitloop in units) {
for (Nperloop in Nperunit) {
sd.fe <- colSds(coeflist.fe[[index]])
sd.re <- colSds(coeflist.re[[index]])
sd.pool <- colSds(coeflist.pool[[index]])
se.fe <- colMeans(selist.fe[[index]])
se.re <- colMeans(selist.re[[index]])
se.pool <- colMeans(selist.pool[[index]])
# opt.fe <- se.fe/sd.fe
# opt.re <- se.re/sd.re
# opt.pool <- se.pool/sd.pool
opt.fe <- sd.fe/se.fe
opt.re <- sd.re/se.re
opt.pool <- sd.pool/se.pool
plot(truecor, opt.fe, pch=19, type="l", col = "#80b1d3", lwd = 3, ylim=c(0, 4), xlim=c(0, 1),
xlab=expression(rho), ylab="Overconfidence",
main=paste(unitloop, "units, ", Nperloop, "observations per unit"))
lines(truecor, opt.re, pch=19, col = "#fb8072", lwd = 3, lty=2)
lines(truecor, opt.pool, pch=19, col = "#bebada", lwd = 3, lty=3)
index <- index + 1
}
}
dev.off()
pdf("CL_Optimism_Sluggish.pdf", width=9, height=8)
par(mfrow=c(length(units), length(Nperunit)), mar=c(4, 4, 2, 0.5), las=1)
index <- 1
for (unitloop in units) {
for (Nperloop in Nperunit) {
sd.fe <- colSds(coeflist.fe[[index]])
sd.re <- colSds(coeflist.re[[index]])
sd.pool <- colSds(coeflist.pool[[index]])
se.fe <- colMeans(selist.fe[[index]])
se.re <- colMeans(selist.re[[index]])
se.pool <- colMeans(selist.pool[[index]])
# opt.fe <- se.fe/sd.fe
# opt.re <- se.re/sd.re
# opt.pool <- se.pool/sd.pool
opt.fe <- sd.fe/se.fe
opt.re <- sd.re/se.re
opt.pool <- sd.pool/se.pool
plot(truecor, opt.fe, pch=19, type="l", col = "#80b1d3", lwd = 3, ylim=c(0, 6), xlim=c(0, 1),
xlab=expression(rho), ylab="Overconfidence",
main=paste(unitloop, "units, ", Nperloop, "observations per unit"))
lines(truecor, opt.re, pch=19, col = "#fb8072", lwd = 3, lty=2)
lines(truecor, opt.pool, pch=19, col = "#bebada", lwd = 3, lty=3)
index <- index + 1
}
}
dev.off()
pdf("CL_Coverage_Sluggish.pdf", width=9, height=8)
par(mfrow=c(length(units), length(Nperunit)), mar=c(4, 4, 2, 0.5), las=1)
index <- 1
for (unitloop in units) {
for (Nperloop in Nperunit) {
coverage.fe <- colMeans(ifelse(lllist.fe[[index]] > b1 | ullist.fe[[index]] < b1, 0, 1))
coverage.re <- colMeans(ifelse(lllist.re[[index]] > b1 | ullist.re[[index]] < b1, 0, 1), na.rm = T)
coverage.pool <- colMeans(ifelse(lllist.pool[[index]] > b1 | ullist.pool[[index]] < b1, 0, 1))
plot(truecor, coverage.fe, pch=19, type="l", col = "#80b1d3", lwd = 3, ylim=c(0, 1), xlim=c(0, 1),
xlab=expression(rho), ylab="Coverage",
main=paste(unitloop, "units, ", Nperloop, "observations per unit"))
lines(truecor, coverage.re, pch=19, col = "#fb8072", lwd = 3, lty=2)
lines(truecor, coverage.pool, pch=19, col = "#bebada", lwd = 3, lty=3)
index <- index + 1
}
}
dev.off()
#quartz.save("CL_Coverage_Sluggish", type = "pdf")
# Plot Optimism
## FE=solid line, RE=dashed line, pooled=dotted line
#quartz(width=9, height=8)
pdf("CL_Optimism_Sluggish.pdf", width=9, height=8)
par(mfrow=c(length(units), length(Nperunit)), mar=c(4, 4, 2, 0.5), las=1)
index <- 1
for (unitloop in units) {
for (Nperloop in Nperunit) {
sd.fe <- colSds(coeflist.fe[[index]])
sd.re <- colSds(coeflist.re[[index]])
sd.pool <- colSds(coeflist.pool[[index]])
se.fe <- colMeans(selist.fe[[index]])
se.re <- colMeans(selist.re[[index]])
se.pool <- colMeans(selist.pool[[index]])
# opt.fe <- se.fe/sd.fe
# opt.re <- se.re/sd.re
# opt.pool <- se.pool/sd.pool
opt.fe <- sd.fe/se.fe
opt.re <- sd.re/se.re
opt.pool <- sd.pool/se.pool
plot(truecor, opt.fe, pch=19, type="l", col = "#80b1d3", lwd = 3, ylim=c(0, 8), xlim=c(0, 1),
xlab=expression(rho), ylab="Overconfidence",
main=paste(unitloop, "units, ", Nperloop, "observations per unit"))
lines(truecor, opt.re, pch=19, col = "#fb8072", lwd = 3, lty=2)
lines(truecor, opt.pool, pch=19, col = "#bebada", lwd = 3, lty=3)
index <- index + 1
}
}
dev.off()
# Plot coverage (how often do our 95% CI's reject beta= truth?):
## FE=solid line, RE=dashed line, pooled=dotted line
#quartz(width=9, height=8)
pdf("CL_Coverage_Sluggish.pdf", width=9, height=8)
par(mfrow=c(length(units), length(Nperunit)), mar=c(4, 4, 2, 0.5), las=1)
index <- 1
for (unitloop in units) {
for (Nperloop in Nperunit) {
coverage.fe <- colMeans(ifelse(lllist.fe[[index]] > b1 | ullist.fe[[index]] < b1, 0, 1))
coverage.re <- colMeans(ifelse(lllist.re[[index]] > b1 | ullist.re[[index]] < b1, 0, 1), na.rm = T)
coverage.pool <- colMeans(ifelse(lllist.pool[[index]] > b1 | ullist.pool[[index]] < b1, 0, 1))
plot(truecor, coverage.fe, pch=19, type="l", col = "#80b1d3", lwd = 3, ylim=c(0, 1), xlim=c(0, 1),
xlab=expression(rho), ylab="Coverage",
main=paste(unitloop, "units, ", Nperloop, "observations per unit"))
lines(truecor, coverage.re, pch=19, col = "#fb8072", lwd = 3, lty=2)
lines(truecor, coverage.pool, pch=19, col = "#bebada", lwd = 3, lty=3)
index <- index + 1
}
}
dev.off()
#quartz.save("CL_Coverage_Sluggish", type = "pdf")
# Plot Optimism
## FE=solid line, RE=dashed line, pooled=dotted line
#quartz(width=9, height=8)
pdf("CL_Optimism_Sluggish.pdf", width=9, height=8)
par(mfrow=c(length(units), length(Nperunit)), mar=c(4, 4, 2, 0.5), las=1)
index <- 1
for (unitloop in units) {
for (Nperloop in Nperunit) {
sd.fe <- colSds(coeflist.fe[[index]])
sd.re <- colSds(coeflist.re[[index]])
sd.pool <- colSds(coeflist.pool[[index]])
se.fe <- colMeans(selist.fe[[index]])
se.re <- colMeans(selist.re[[index]])
se.pool <- colMeans(selist.pool[[index]])
# opt.fe <- se.fe/sd.fe
# opt.re <- se.re/sd.re
# opt.pool <- se.pool/sd.pool
opt.fe <- sd.fe/se.fe
opt.re <- sd.re/se.re
opt.pool <- sd.pool/se.pool
plot(truecor, opt.fe, pch=19, type="l", col = "#80b1d3", lwd = 3, ylim=c(0, 10), xlim=c(0, 1),
xlab=expression(rho), ylab="Overconfidence",
main=paste(unitloop, "units, ", Nperloop, "observations per unit"))
lines(truecor, opt.re, pch=19, col = "#fb8072", lwd = 3, lty=2)
lines(truecor, opt.pool, pch=19, col = "#bebada", lwd = 3, lty=3)
index <- index + 1
}
}
dev.off()
# Plot coverage (how often do our 95% CI's reject beta= truth?):
## FE=solid line, RE=dashed line, pooled=dotted line
#quartz(width=9, height=8)
pdf("CL_Coverage_Sluggish.pdf", width=9, height=8)
par(mfrow=c(length(units), length(Nperunit)), mar=c(4, 4, 2, 0.5), las=1)
index <- 1
for (unitloop in units) {
for (Nperloop in Nperunit) {
coverage.fe <- colMeans(ifelse(lllist.fe[[index]] > b1 | ullist.fe[[index]] < b1, 0, 1))
coverage.re <- colMeans(ifelse(lllist.re[[index]] > b1 | ullist.re[[index]] < b1, 0, 1), na.rm = T)
coverage.pool <- colMeans(ifelse(lllist.pool[[index]] > b1 | ullist.pool[[index]] < b1, 0, 1))
plot(truecor, coverage.fe, pch=19, type="l", col = "#80b1d3", lwd = 3, ylim=c(0, 1), xlim=c(0, 1),
xlab=expression(rho), ylab="Coverage",
main=paste(unitloop, "units, ", Nperloop, "observations per unit"))
lines(truecor, coverage.re, pch=19, col = "#fb8072", lwd = 3, lty=2)
lines(truecor, coverage.pool, pch=19, col = "#bebada", lwd = 3, lty=3)
index <- index + 1
}
}
dev.off()
#quartz.save("CL_Coverage_Sluggish", type = "pdf")
# Plot Optimism
## FE=solid line, RE=dashed line, pooled=dotted line
#quartz(width=9, height=8)
pdf("CL_Optimism_Sluggish.pdf", width=9, height=8)
par(mfrow=c(length(units), length(Nperunit)), mar=c(4, 4, 2, 0.5), las=1)
index <- 1
for (unitloop in units) {
for (Nperloop in Nperunit) {
sd.fe <- colSds(coeflist.fe[[index]])
sd.re <- colSds(coeflist.re[[index]])
sd.pool <- colSds(coeflist.pool[[index]])
se.fe <- colMeans(selist.fe[[index]])
se.re <- colMeans(selist.re[[index]])
se.pool <- colMeans(selist.pool[[index]])
# opt.fe <- se.fe/sd.fe
# opt.re <- se.re/sd.re
# opt.pool <- se.pool/sd.pool
opt.fe <- sd.fe/se.fe
opt.re <- sd.re/se.re
opt.pool <- sd.pool/se.pool
plot(truecor, opt.fe, pch=19, type="l", col = "#80b1d3", lwd = 3, ylim=c(0, 10), xlim=c(0, 1),
xlab=expression(rho), ylab="Overconfidence",
main=paste(unitloop, "units, ", Nperloop, "observations per unit"))
lines(truecor, opt.re, pch=19, col = "#fb8072", lwd = 3, lty=2)
lines(truecor, opt.pool, pch=19, col = "#bebada", lwd = 3, lty=3)
index <- index + 1
}
}
dev.off()
pdf("CL_Optimism_Sluggish.pdf", width=9, height=8)
par(mfrow=c(length(units), length(Nperunit)), mar=c(4, 4, 2, 0.5), las=1)
index <- 1
for (unitloop in units) {
for (Nperloop in Nperunit) {
sd.fe <- colSds(coeflist.fe[[index]])
sd.re <- colSds(coeflist.re[[index]])
sd.pool <- colSds(coeflist.pool[[index]])
se.fe <- colMeans(selist.fe[[index]])
se.re <- colMeans(selist.re[[index]])
se.pool <- colMeans(selist.pool[[index]])
# opt.fe <- se.fe/sd.fe
# opt.re <- se.re/sd.re
# opt.pool <- se.pool/sd.pool
opt.fe <- sd.fe/se.fe
opt.re <- sd.re/se.re
opt.pool <- sd.pool/se.pool
plot(truecor, opt.fe, pch=19, type="l", col = "#80b1d3", lwd = 3, ylim=c(0, 6), xlim=c(0, 1),
xlab=expression(rho), ylab="Overconfidence",
main=paste(unitloop, "units, ", Nperloop, "observations per unit"))
lines(truecor, opt.re, pch=19, col = "#fb8072", lwd = 3, lty=2)
lines(truecor, opt.pool, pch=19, col = "#bebada", lwd = 3, lty=3)
index <- index + 1
}
}
dev.off()
load("/Users/alikagalwala/Dropbox/CSTS-Panel Papers/JOP RRR/Replication Files/Clark and Linzer 2015/Sluggish/2000sims_sluggish.RData")
normal.pool <- matrix(data = NA, nrow = 99, ncol = 6)
colnames(normal.pool) <- c("units", "within_units", "corr", "p_value", "normal_90", "normal_95")
normal.fe <- matrix(data = NA, nrow = 99, ncol = 6)
colnames(normal.fe) <- c("units", "within_units", "corr", "p_value", "normal_90", "normal_95")
normal.re <- matrix(data = NA, nrow = 99, ncol = 6)
colnames(normal.re) <- c("units", "within_units", "corr", "p_value", "normal_90", "normal_95")
index <- 1
j <- 1
for (unitloop in units) {
for (Nperloop in Nperunit) {
for (i in 1:length(truecor)) {
# pooled
normal.pool[index, 1] <- unitloop
normal.pool[index, 2] <- Nperloop
normal.pool[index, 3] <- truecor[i]
normal.pool[index, 4] <- shapiro.test(coeflist.pool[[j]][ , i])$p.value
normal.pool[index, 5] <- ifelse(shapiro.test(coeflist.pool[[j]][ , i])$p.value > 0.1, "normal", "not normal")
normal.pool[index, 6] <- ifelse(shapiro.test(coeflist.pool[[j]][ , i])$p.value > 0.05, "normal", "not normal")
# fe
normal.fe[index, 1] <- unitloop
normal.fe[index, 2] <- Nperloop
normal.fe[index, 3] <- truecor[i]
normal.fe[index, 4] <- shapiro.test(coeflist.fe[[j]][ , i])$p.value
normal.fe[index, 5] <- ifelse(shapiro.test(coeflist.fe[[j]][ , i])$p.value > 0.1, "normal", "not normal")
normal.fe[index, 6] <- ifelse(shapiro.test(coeflist.fe[[j]][ , i])$p.value > 0.05, "normal", "not normal")
# re
normal.re[index, 1] <- unitloop
normal.re[index, 2] <- Nperloop
normal.re[index, 3] <- truecor[i]
normal.re[index, 4] <- shapiro.test(coeflist.re[[j]][ , i])$p.value
normal.re[index, 5] <- ifelse(shapiro.test(coeflist.re[[j]][ , i])$p.value > 0.1, "normal", "not normal")
normal.re[index, 6] <- ifelse(shapiro.test(coeflist.re[[j]][ , i])$p.value > 0.05, "normal", "not normal")
index <- index + 1
print(j)
}
j <- j + 1
}
}
normal.pool
normal.fe
normal.re
